/// \file solvesystems.h
/// \brief This file contains many general functions for matrix algebra
///
///
///
/// \author  Wim van Ackooij
/// \date  may 2007
///
#ifndef __SOLVESYSTEMS_HPP__
#define __SOLVESYSTEMS_HPP__

#ifndef HAVE_MDMALLOC
#ifndef GN_MALLOC
#define GN_MALLOC(tab, type, size, routine){						\
  tab = (type*)malloc(size*sizeof(type));						\
  if (tab==NULL) { printf("\n\n !!! ERREUR !!! [%s] : erreur d'allocation memoire\n\n",routine); exit(1); }	\
}
#endif

#ifndef GN_CALLOC
#define GN_CALLOC(tab, type, size, routine){												\
  tab = (type*)calloc(size,sizeof(type));												\
  if (tab==NULL) { printf("\n\n !!! ERREUR !!! [%s] : erreur d'allocation memoire\n\n",routine); exit(1); }	\
}
#endif
#endif

#ifdef HAVE_MDMALLOC
extern "C" {
#include "projects/libv0/src/netcdfStruct/mdmalloc.h"
}
#endif

#include <assert.h>
#include <iostream>
#include <sstream>
#include <cmath>
#include <vector>
#include <iomanip>
#include <string>

#include <blitz/array.h>
#include <boost/date_time/gregorian/gregorian.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>

#include "src/common/general.h"

// Functions

/// \brief  Invert the matrix using the Gauss Algorithm
///		N.B. : Gauss algorithm is known for its poor numerical stability
/// \param A	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \return C	(double**)  the inverse of A
double** gaussinv(double**, int); // Inversion de Matrice avec Algo de Gauss

/// \brief  Multiply square matrix with a vector
/// \param A	(double**)	the Matrix
/// \param x	(double*)	the vector
/// \param nline	(int)		the number of lines / columns
/// \return y	(double*)   the vector y = Ax
double* matvect(double**&, double*&, int);    // calcule A b

/// \brief  Multiply transposed (m times n) matrix with a vector
/// \param A	(double**)	the Matrix
/// \param x	(double*)	the vector (of size m)
/// \param m	(int)		the number of lines
/// \param n	(int)		the number of columns
/// \return y	(double*)   the vector y = A^Tx
double* matTvect(double**&, double*&, int, int);   // calcule A^T b A est m x n, b est m x 1

/// \brief  Multiply square matrices
/// \param A	(double**)	the Matrix
/// \param B	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \return C	(double**)   the vector C = A*B
double** matmat(double**&, double**&, int n);   // calcule AB

/// \brief  Multiply transposed matrix A with matrix B (both of size m times n)
/// \param A	(double**)	the Matrix
/// \param B	(double**)	the Matrix
/// \param m	(int)		the number of lines
/// \param n	(int)		the number of columns
/// \return C	(double**)   the vector C = A^T*B (size n times n)
double** matTmat(double**&, double**&, int m, int n);   // calcule A^T B, les deux sont matrices m x n

/// \brief  Multiply square matrix B before and after with A
/// \param A	(double**)	the Matrix
/// \param B	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \return C	(double**)   the vector C = A*B*A
double** matprepost(double**&, double**&, int);   // calcul ABA

/// \brief  Multiply square matrix A with a diagonal D and then with square matrix B
/// \param A	(double**)	the Matrix
/// \param D	(double*)	the diagonal matrix
/// \param B	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \return C	(double**)   the vector C = A*D*B
double** matdiag(double**&, double*&, double**&, int);    // calcul ADB, where D is a diagonal matrix

/// \brief  Multiply square matrix A with a diagonal D and then with transposed square matrix B
/// \param A	(double**)	the Matrix
/// \param D	(double*)	the diagonal matrix
/// \param B	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \return C	(double**)   the vector C = A*D*B^T
double** matdiagtranspose(double**&, double*&, double**&, int);    // calcul ADB^T, where D is a diagonal matrix

/// \brief  Solve the matrix vector system Ax = b
///		Uses the Richardson-Tanabe variant of the algorithm found in :
///		R. Zorgati, W. van Ackooij and M. Lambert.
///		Stochastic Matrice and Lp norms : New Algorithms for solving Ill-conditioned linear system of Equations.
///		European Series in Applied and Industrial Mathematics (ESAIM) proc vol 18 ; 70-86. 2007.
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* zrg_rt(double**&, double*&, int);

/// \brief  Solve the matrix vector system Ax = b
///		Uses the Schultz-Hotelling-Bodewig variant of the algorithm found in :
///		R. Zorgati, W. van Ackooij and M. Lambert.
///		Stochastic Matrice and Lp norms : New Algorithms for solving Ill-conditioned linear system of Equations.
///		European Series in Applied and Industrial Mathematics (ESAIM) proc vol 18 ; 70-86. 2007.
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* zrg_shb(double**&, double*&, int);

/// \brief  The fortran 77 sign function
/// \param  x 	(double) 	some double
/// \param  y 	(double) 	some double
/// \return z	(double)	returns value of x, with sign of y, i.e. z = abs(x)sign(y)
double sign(double, double);

/// \brief  The fortran 77 pythag function
/// \param  x 	(double) 	some double
/// \param  y 	(double) 	some double
/// \return z	(double)	returns z = sqrt(x^2 + y^2)
double pythag(double,double);

/// \brief QR decomposition
///	     Constructs the QR decomposition of a[1..n][1..n].
///	     The upper triangular matrix R is returned in the upper triangle of a, except for the diagonal elements of R which are returned in d[1..n].
///	     The orthogonal matrix Q is represented as a product of n Householder matrices Q_1  x ... x Q_n . 1 , where Q_j =1 - u_j tensor u_j / c_j.
///	     The ith component of u_j is zero for i = 1 ... j - 1
///	     while the nonzero components are returned in a[i][j] for i = j ... n.
///	     sing returns as true (1) if singularity is encountered during the decomposition, but the decomposition is still completed in this case; otherwise it returns false (0).
/// \param A	(double**)	the Matrix to decompose
/// \param n	(int)		the size of A
/// \param R	(double**)	the Matrix R from A = QR
/// \param c	(double*)	element to build Q
/// \param d	(double*)	the diagonal of matrix R
/// \param is_sing(bool)	did we encounter singularities
void QRdcmp(double** const, int n, double**&, double*&, double*&, bool);

/// \brief  Solve the matrix vector system Ax = b using QR decomposition
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* QRsolve(double** const, double* const, int);

/// \brief  Solve the matrix vector system Ax = b using QR decomposition
/// 		We should call QRdcmp first and use those elements as input here
/// \param b	(double*)	the right member
/// \param n	(int)		the size of A
/// \param R	(double**)	the Matrix R from A = QR
/// \param c	(double*)	element to build Q
/// \param d	(double*)	the diagonal of matrix R
/// \param is_sing(bool)	did we encounter singularities
double* QRsolve(double*, int, double** const, double* const, double* const, bool); //Should run QRdcmp first on the matrix !!!

/// \brief  Solve the matrix vector system Ax = b using QR decomposition
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \param x	(double*)   the vector x, s.t. Ax = b, This should be already allocated !!!
void QRsolve(double**&, double*&, int, double*&);

/// \brief LU decomposition
/// 		Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise
/// 		permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above;
/// 		indx[1..n] is an output vector that records the row permutation effected by the partial
/// 		pivoting; d is output as +/-1 depending on whether the number of row interchanges was even
/// 		or odd, respectively. This routine is used in combination with lubksb to solve linear equations
/// 		or invert a matrix.
/// \param A	(double**)	the Matrix to decompose
/// \param U	(double**)	the Matrix such that A = LU
/// \param n	(int)		the size of A
/// \param indx	(int*)	records the row permutation effected by the partial pivoting
/// \param d	(int)		+1 if number of row changes was even, -1 if this was odd
void LUdcmp(double**&, double**&, int n, int*&, int);

/// \brief  Solve the matrix vector system Ax = b using LU decomposition
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* LUsolve(double**&, double*&, int);

/// \brief SVD decomposition
/// 		Given a matrix a[1..m][1..n], this routine computes its singular value
/// 		decomposition, A = U.W.V^T.  The matrix U replaces a on output.  The diagonal
/// 		matrix of singular values W is output as a vector w[1..n].  The matrix V (not
/// 		the transpose VT) is output as v[1..n][1..n].
/// 		Given a matrix A m x n, computes the SVD, A = U W V^T,
/// 		where W is output as a vector, as it is a diagonal matrix
/// \param A	(double**)	the Matrix
/// \param m	(int)		the number of lines
/// \param n	(int)		the number of columns
/// \param U	(double**)	the Matrix U
/// \param w	(double*)	the diagonal matrix W
/// \param V	(double**)	the Matrix V
void SVDdcmp(double**&, int,int, double**&, double*&, double**&);

/// \brief  Solve the matrix vector system Ax = b using SVD decomposition
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param m	(int)		the number of lines
/// \param n	(int)		the number of columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* SVDsolve(double**&, double*&, int, int);

/// \brief SVD decomposition
///		blitz version of SVDdcmp
/// \author X Warin
/// \date jan 2001 , aout 2002, janvier 2003, mars 2006 (version C++ indicage 0)
/// 		Given a matrix a[1..m][1..n], this routine computes its singular value
/// 		decomposition, A = U.W.V^T.  The matrix U replaces a on output.  The diagonal
/// 		matrix of singular values W is output as a vector w[1..n].  The matrix V (not
/// 		the transpose VT) is output as v[1..n][1..n].
/// 		Given a matrix A m x n, computes the SVD, A = U W V^T,
/// 		where W is output as a vector, as it is a diagonal matrix
///	Attention : This version destroys A
/// \param A	(blitz::Array<double,2>)	the Matrix A (as output this is U)
/// \param w	(blitz::Array<double,1>)	the diagonal matrix W
/// \param V	(blitz::Array<double,2>)	the Matrix V
void svdcmp( blitz::Array<double, 2> & a , blitz::Array<double,1> & w , blitz::Array<double,2> &v);

/// \brief  Solve the matrix vector system Ax = b using SVD decomposition
///		Input should come from the svddcmp function
/// \author X Warin
/// \date 08/01/2002, mars 2006
/// \param U	(blitz::Array<double,2>)	the Matrix U
/// \param w	(blitz::Array<double,1>)	the diagonal matrix W
/// \param V	(blitz::Array<double,2>)	the Matrix V
/// \param b	(blitz::Array<double,1>)	the right member
/// \param x	(blitz::Array<double,1>)   	the vector x, s.t. Ax = b
void svdksb(blitz::Array<double,2> & u, blitz::Array< double,1> & w, blitz::Array<double,2> & v, blitz::Array<double,1> & b, blitz::Array<double,1> & x);

/// \brief 	Checks if a matrix is a diagonal matrix
/// \param A	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \return iD	(bool)	true / false
bool is_diagMat(double**&, int);

/// \brief  Solve the matrix vector system Ax = b assuming that A is diagonal
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* Diagsolve(double** const, double* const, int);

/// \brief Choleski decomposition : A = CC^T
/// \param A	(double**)	the Matrix
/// \param nline	(int)		the number of lines / columns
/// \param isOk	(bool)	is the matrix positive semi definite ? (true / false)
/// \return C	(double**)  the matrix such that A = C C^T
double** Choldcmp( double**, int, bool&);

/// \brief  Solve the matrix vector system Ax = b using choleski decomposition
///		Attention : this version destroys A
/// \param A	(double**)	the Matrix
/// \param b	(double*)	the right member
/// \param nline	(int)		the number of lines / columns
/// \return x	(double*)   the vector x, s.t. Ax = b
double* Cholsolve( double**& , double*&, int );   //Attention detruit A

/// \brief Choleski decomposition : A = CC^T
///		blitz version of choldcmp
/// \param A	(blitz::Array<double,2>)	the Matrix
/// \param isOk	(bool)	is the matrix positive semi definite ? (true / false)
/// \return C	(blitz::Array<double,2>)  the matrix such that A = C C^T
blitz::Array<double,2> Choldcmp( blitz::Array<double,2> & A, bool& isOk);


#endif
